% By Martin Andreasen, 2009
% This function computes second order analytical derivatives 
% of the pricing kernel. The integer "R_function" has the following properties:
% 1  - if log-transformation of bond prices is used and 
% 0  - if no log-transformation of bond prices is used

function [M,Mx,Mxp,My,Myp,...
    Mxx,Mxxp,Mxy,Mxyp,Mxpx,Mxpxp,Mxpy,Mxpyp,Myx,Myxp,Myy,Myyp,Mypx,Mypxp,Mypy,Mypyp] = ...
    Anal_PricingKernel_derivatives(x,xp,y,yp,order_app,R_function);

% Some integers
nx    = size(x,2);
ny    = size(y,2);

%Define parameters
syms BETTA B CHI TAU THETA DELTA ALFA PHI PHIzero ZI ETA RHOR PHIpai PHIy  ...
     A_aa     A_az     A_ad     A_adelta...
     A_za     A_zz     A_zd     A_zdelta...
     A_da     A_dz     A_dd     A_ddelta...
     A_deltaa A_deltaz A_deltad A_deltadelta ...
     Kss Rss OUTPUTss PAIss MUZss AA;

%Define variables 
syms r_cu   c_cu   pai_cu   w_cu    l_cu    b_cu      evf_cu   vf_cu...
     r_cup  c_cup  pai_cup  w_cup   l_cup   b_cup     evf_cup  vf_cup...
     r_ba1  c_ba1  a_cu     muz_cu  d_cu    delta_cu  epsR_cu  ...
     r_ba1p c_ba1p a_cup    muz_cup d_cup   delta_cup epsR_cup ;

% Some static relations which we substitute into f
% EQ 4: Firms FOC for labour
%mc_cu  = w_cu /((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);
%mc_cup = w_cup/((1-THETA)*a_cup*Kss^THETA*l_cup^-THETA);

% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(muz_cup^CHI*vf_cup))^ALFA...   
    *((((1-exp(-TAU*(c_cup-B*c_cu/muz_cup)))/TAU)^(CHI-1))*exp(-TAU*(c_cup-B*c_cu/muz_cup)))...
    /((((1-exp(-TAU*(c_cu-B*c_ba1/muz_cu)))/TAU)^(CHI-1))*exp(-TAU*(c_cu-B*c_ba1/muz_cu)))*muz_cup^(CHI-1);  

% EQ 7: Output
%output_cu  = (c_cu  + delta_cu *Kss)/(1-ZI/2*(pai_cu /PAIss-1)^2);
%output_cup = (c_cup + delta_cup*Kss)/(1-ZI/2*(pai_cup/PAIss-1)^2);


% The pricing kernel
M = mReal_cup/pai_cup;

%Make f a function of the logarithm of the state and control vector
if R_function == 1
    y_log  = [r_cu   c_cu   pai_cu   w_cu    l_cu    b_cu      evf_cu   vf_cu];
    yp_log = [r_cup  c_cup  pai_cup  w_cup   l_cup   b_cup     evf_cup  vf_cup];
    x_log  = [r_ba1  c_ba1  a_cu     muz_cu  d_cu    delta_cu  ];
    xp_log = [r_ba1p c_ba1p a_cup    muz_cup d_cup   delta_cup ];

    M = subs(M, [x_log,y_log,xp_log,yp_log], exp([x_log,y_log,xp_log,yp_log]));
end
% For Extended Path
M_ep = char(M);
M_ep = strrep(M_ep,'*','.*');
M_ep = strrep(M_ep,'/','./');
M_ep = strrep(M_ep,'^','.^')

% The first order derivatives of the pricing kernel
% Recall we only need the first order derivatives of the pricing kernel to get
% bond prices up to second order
if order_app > 1
    Mx  = jacobian(M,x);
    Mxp = jacobian(M,xp);
    My  = jacobian(M,y);
    Myp = jacobian(M,yp);
else
    Mx  = zeros(1,nx);
    Mxp = zeros(1,nx);
    My  = zeros(1,ny);
    Myp = zeros(1,ny);
end

% The second order derivatives of the pricing kernel
% Recall we only need the second order derivatives of the pricing kernel to get
% bond prices up to third order
if order_app > 2 
    % With respect to Mx
    Mxx   = reshape(jacobian(Mx(:) ,x) ,nx,nx);
    Mxxp  = reshape(jacobian(Mx(:) ,xp),nx,nx);
    Mxy   = reshape(jacobian(Mx(:) ,y) ,nx,ny);
    Mxyp  = reshape(jacobian(Mx(:) ,yp),nx,ny);

    % With respect to Mxp
    Mxpx  = reshape(jacobian(Mxp(:),x) ,nx,nx);
    Mxpxp = reshape(jacobian(Mxp(:),xp),nx,nx);
    Mxpy  = reshape(jacobian(Mxp(:),y) ,nx,ny);
    Mxpyp = reshape(jacobian(Mxp(:),yp),nx,ny);

    % With respect to My
    Myx   = reshape(jacobian(My(:) ,x) ,ny,nx);
    Myxp  = reshape(jacobian(My(:) ,xp),ny,nx);
    Myy   = reshape(jacobian(My(:) ,y) ,ny,ny);
    Myyp  = reshape(jacobian(My(:) ,yp),ny,ny);

    % With respect to Myp
    Mypx  = reshape(jacobian(Myp(:),x ),ny,nx);    
    Mypxp = reshape(jacobian(Myp(:),xp),ny,nx);
    Mypy  = reshape(jacobian(Myp(:),y),ny,ny);
    Mypyp = reshape(jacobian(Myp(:),yp),ny,ny);
else
    % With respect to Mx
    Mxx   = zeros(nx,nx);
    Mxxp  = zeros(nx,nx);
    Mxy   = zeros(nx,ny);
    Mxyp  = zeros(nx,ny);

    % With respect to Mxp
    Mxpx  = zeros(nx,nx);
    Mxpxp = zeros(nx,nx);
    Mxpy  = zeros(nx,ny);
    Mxpyp = zeros(nx,ny);

    % With respect to My
    Myx   = zeros(ny,nx);
    Myxp  = zeros(ny,nx);
    Myy   = zeros(ny,ny);
    Myyp  = zeros(ny,ny);

    % With respect to Myp
    Mypx  = zeros(ny,nx);    
    Mypxp = zeros(ny,nx);    
    Mypy  = zeros(ny,ny);    
    Mypyp = zeros(ny,ny);
end
